Introducción al manejo de datos territoriales
La librería sf en R implementa un manejo eficiente y estandarizado de datos espaciales utilizando el modelo de geometrías de características simples (o simple features), un estándar internacional definido por ISO 19125 y desarrollado por el Open Geospatial Consortium (OGC). Este modelo es ampliamente utilizado en herramientas de software geográfico como GeoJSON, ArcGIS, QGIS, PostGIS, MySQL Spatial Extensions, y Microsoft SQL Server, facilitando la interoperabilidad y el análisis de datos geoespaciales en diversos entornos.
En R, el paquete sf proporciona una representación poderosa y flexible de datos espaciales vectoriales. Los objetos principales que maneja son extensiones de los clásicos data frames o tibbles, y están diseñados para facilitar la integración con otras herramientas de análisis de datos en R. Cada objeto del tipo sf incluye una columna especial denominada simple feature geometry list column, que almacena la geometría de cada observación en un formato tipo lista. Esta geometría se asocia con una o más variables adicionales (atributos) que describen cada entidad geográfica, lo que convierte cada fila en una simple feature o característica simple, compuesta tanto por sus atributos como por su estructura geométrica.
Además, sf introduce una serie de funciones que permiten realizar transformaciones geométricas, operaciones topológicas y análisis espaciales de manera eficiente, aprovechando bibliotecas subyacentes como GEOS, GDAL y PROJ. Este enfoque no solo facilita la manipulación de datos geoespaciales, sino que también asegura una alta compatibilidad con los formatos y herramientas geográficas más utilizados a nivel global.
Para comenzar, instalamos la librería y la cargamos a nuestro entorno de trabajo:
install.packages("sf")
library(sf)Tipos de datos espaciales
Es importante distinguir entre los dos tipos principales de datos espaciales: datos vectoriales y datos raster. Ambos tipos de datos permiten representar información geográfica, pero lo hacen de manera muy diferente.
Datos vectoriales
Los datos vectoriales son el enfoque utilizado por la librería sf y se basan en representaciones geométricas precisas de objetos espaciales a través de puntos, líneas y polígonos. Estos son comúnmente llamados features o características simples. Los tipos principales de geometrías que manejan los datos vectoriales son:
Puntos (Point)
Representan ubicaciones específicas en un espacio de coordenadas. Cada punto está definido por un par de coordenadas (x, y) en el espacio 2D, o (x, y, z) en el espacio 3D.
Para comprender mejor la lógica, vamos a crear un punto en el espacio bidimensional de las variables x e y. En primer lugar, definimos un vector con las coordenadas que asignaremos a este punto:
coordenadas = c(5,3)
coordenadas## [1] 5 3Es decir, el punto estará ubicado en el plano carteriano en una localización definida por x = 5 e y = 3.
A continuación, utilizamos la función st_point para convertir a este vector en un objeto espacial, a partir de esas coordinadas definidas previamente:
punto <- st_point(coordenadas)
punto## POINT (5 3)Si graficamos este objeto, observamos lo siguiente:
plot(punto, axes = TRUE)
Si bien este punto está localizado en el plano carteriano a partir de los valores 5 (para x) y 3 (para y), la lógica aplica cuando estos valores no representan ejes abstractos sino un sistema de coordenadas dado por la latitud y la longitud. En estos casos, como veremos más adelante, el punto estaría localizado sobre el espacio geográfico.
También podríamos contar con una colección de puntos dentro de un mismo objeto. Para ver esto, en primer lugar definimos un vector de valores para x junto con un vector de valores para y, haciendo:
px = c(5, 4, 2)
py = c(3, 5, 4)Ahora, combinamos este conjunto de valores en un mismo data frame mediante la función cbind():
coordenadas <- cbind(px, py)
coordenadas## px py
## [1,] 5 3
## [2,] 4 5
## [3,] 2 4Como se puede apreciar, el primer punto estará ubicado en la localización definida por los valores x = 5 e y = 3. El segundo en la localización dada por x = 4 e y = 5, y el tercer punto en la localización definida por el par ordenado x = 2 e y = 4.
Al igual que en el caso de un solo punto, convertimos a este data frame en un objeto espacial, pero ahora utilizando la función st_multipoint():
puntos <- st_multipoint(coordenadas)
puntos## MULTIPOINT ((5 3), (4 5), (2 4))Graficando este objeto:
plot(puntos, axes = TRUE)
Líneas
Las líneas son secuencias de puntos conectados que representan entidades lineales como caminos, ríos o vías de tren. Para comprender la lógica de este tipo de datos vectoriales, comencemos por definir un conjunto de coordenadas para la línea “l1”.
l1x <- c(0, 3, 5, 8, 10)
l1y <- c(0, 4, 4, 8, 10)
coordenadas_l1 <- cbind(l1x, l1y)
coordenadas_l1## l1x l1y
## [1,] 0 0
## [2,] 3 4
## [3,] 5 4
## [4,] 8 8
## [5,] 10 10Es decir que la línea l1 partirá el origen (en donde x = 0 e y = 0), luego se dirigirá hacia el punto dado por x = 3 e y = 4, luego pasará por x = 5 e y = 4, y proseguirá su viaje hasta arribar a destino en la localización dada por el punto x = 10 e y = 10.
A partir de estas coordenadas, vamos a crear un objeto espacial utilizando la función st_linestring():
l1 <- st_linestring(coordenadas_l1)
plot(l1, axes = TRUE)
A continuación, podemos repetir el mismo procedimiento para crear una línea diferente en el mismo espacio cartesiano:
l2x <- c(2, 3, 12)
l2y <- c(4, 8, 5)
coordenadas_l2 <- cbind(l2x, l2y)
l2 <- st_linestring(coordenadas_l2)
plot(l2, axes = TRUE)
Finalmente, podemos combinar ambas líneas en un mismo objeto espacial, haciendo uso de la función st_multilinestring():
lineas <- st_multilinestring(list(l1, l2))
plot(lineas, axes = TRUE)
Polígonos
Los polígonos son representaciones formadas por líneas que encierran un área, como límites de países, lagos o áreas protegidas. Para comprender la lógica de este tipo de datos vectoriales, comencemos por definir un conjunto de coordenadas para un polígono imaginario en el plano cartesiano.
x <- c(5, 10, 10, 6, 5)
y <- c(1, 2, 5, 4, 1)
coordenadas <- cbind(px, py)
coordenadas## px py
## [1,] 5 3
## [2,] 4 5
## [3,] 2 4La conformación del polígono parte de la lectura de esta lista de localizaciones en sentido antihorario. Es decir, la línea que define el contorno del polígono partirá, por ejemplo, del punto dado por x = 5 e y = 1. De allí se dirigirá hacia el punto localizado en x = 10 e y = 2. La línea continuará su viaje por cada uno de estos puntos que marcan los vértices del polígono, para arribar finalmente, de nuevo al punto de partida, dado por x = 5 e y = 1. Para convertir esta lista de coordenadas en un objeto espacial, utilizamos la función st_polygon():
poligono <- st_polygon(list(coordenadas))Si graficamos este objeto, tenemos que:
plot(poligono, axes = T)
Colecciones geométricas
Combinaciones de puntos, líneas y/o polígonos en un solo objeto. Cada una de estas geometrías puede estar acompañada de atributos, como por ejemplo el nombre de una ciudad para un punto, la longitud de una ruta para una línea, o la población de un área para un polígono. Estos atributos permiten realizar análisis espaciales avanzados y crear mapas temáticos. Para armar un objeto espacial que consista en una colección de los puntos, líneas y polígonos construidos anteriormente, utilizamos la función st_geometrycollection().
coleccion <- st_geometrycollection(list(puntos, lineas, poligono))
plot(st_sfc(coleccion), axes = TRUE)
Datos raster
A diferencia de los datos vectoriales, los datos raster dividen el espacio geográfico en una cuadrícula regular de celdas o píxeles. Cada celda contiene un valor numérico que representa información continua o discreta, como la elevación de un terreno, la temperatura en una región o la clasificación del uso del suelo. Los datos raster son ideales para representar fenómenos que varían de manera continua en el espacio, como el clima, la geología o las imágenes satelitales.
A pesar de que la librería sf está centrada en datos vectoriales, es común trabajar en conjunto con datos raster, utilizando otras librerías en R como raster o terra para gestionar y analizar este tipo de datos. Las herramientas modernas permiten realizar análisis espaciales combinados, aprovechando tanto los datos vectoriales (por ejemplo, límites políticos o redes de transporte) como los datos raster (por ejemplo, modelos de elevación digital o imágenes satelitales).
Estructuralmente, los rásters consisten en matrices de valores. Para comprender mejor la lógica de este tipo de objetivos, en primer lugar, construimos una matriz de la siguiente manera:
matriz <- matrix(c(1, 2, 3, 4, 2, NA, 2, 2, 3, 3, 3, 1), ncol = 4, nrow = 3, byrow = TRUE)
matriz## [,1] [,2] [,3] [,4]
## [1,] 1 2 3 4
## [2,] 2 NA 2 2
## [3,] 3 3 3 1A continuación, vamos a utilizar la librería raster para convertir esta matriz en un objeto espacial. Primero, instalamos la librería:
install.packages("raster")Posteriormente la cargamos al entorno de trabajo, y utilizamos la función raster() para lograr el objetivo:
library(raster)## Loading required package: spr <- raster(matriz)
plot(r)
Como se observa en el gráfico, el raster contiene diferentes valores para cada celda. Los valores más bajos se representan en color gris, en tanto que lo valores más altos se representan en color verde. Las celdas representan la estructura de la matriz, en donde, por ejemplo, en la tercera fila (la de abajo) y la tercera columna (la de la derecha) el valor de la celda es igual a 1, mientras que en la primera fila (la de arriba) y la tercera columna el valor asociado es igual a 4. La intersección entre la segunda fila y la segunda columna de la matriz tiene asociado el valor NA, que se representa como una celda en blanco en el gráfico. En R, NA representa un valor faltante o no disponible (Not Available). Es la forma en que R indica que no tiene un valor para una observación en particular, lo que puede suceder por diversas razones, como datos incompletos, errores de entrada o valores omitidos.
Extrayendo datos geográficos abiertos
La comunidad de R ha generado diferentes herramientas para la extracción de datos libres y abiertos de naturaleza geográfica. Una excelente fuente de información de este tipo es Openstreetmap (OSM), que consiste en un proyecto colaborativo que tiene como objetivo crear un mapa mundial libre y editable por cualquier persona. Se basa en la filosofía de crowdsourcing, lo que significa que los datos del mapa son generados y mantenidos por voluntarios de todo el mundo, quienes contribuyen con información sobre calles, edificios, parques, ríos, y muchos otros elementos geográficos. Para compensar la utilización gratuita de información de la cual gozaremos gracias al esfuerzo desinteresado de miles de personas que aportan información al mapa, recomendamos que las personas se creen un usuario en OSM y naveguen por el mapa, contribuyendo a sumar información de sus entornos y a elevar la calidad del mapa colaborativo. Esto puede realizarse en el www.openstreetmap.org/. El siguiente link contiene una guia paso a paso para comenzar a colaborar con el mapeo.
La librería osmdata permite extraer datos desde la API de OSM, por lo tanto, comenzamos por instalarla. También instalaremos la librería nominatimlite, que necesitaremos para definir entorno geográficos conocidos sólo a partir de su nombre:
install.packages("osmdata")
install.packages("nominatimlite")A continuación, vamos a definir un entorno geográfico (un polígono) sobre el cuál se acotará la extracción de información de OSM. Este polígono vendrá definido por los límites administrativos de la Municipalidad de Córdoba:
library(nominatimlite)
library(osmdata)
area_de_estudio <- geo_lite_sf(address = "Cordoba, Argentina", points_only = F)Veamos un poco más de cerca el polígono area_de_estudio:
area_de_estudio## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -64.3101 ymin: -31.52513 xmax: -64.05725 ymax: -31.30754
## Geodetic CRS: WGS 84
## # A tibble: 1 × 3
## query address geometry
## * <chr> <chr> <POLYGON [°]>
## 1 Cordoba, Argentina Córdoba, Municipio de Córdoba, P… ((-64.3101 -31.36306, -6…Ya entraremos en más detalles respecto a este objeto, pero por lo pronto, notemos que tiene la estructura de un data frame, con una sola fila en donde se observan las siguientes columnas: query, address, geometry. La última columna es la que nos interesa en este momento, dado que es una lista de coordenadas geográficas que nos permitirá visualizar este objeto en un mapa.
Para ello, vamos a hacer uso de la librería mapview, que permite la construcción de mapas interactivos de una manera muy sencilla. Instalamos esta librería:
install.packages("mapview")Y a continuación, construimos el mapa:
library(mapview)
mapview(area_de_estudio)Al hacer click en el polígono notamos que se abre una ventana con el resto de la información contenida en las columnas restantes. Esta es la principal ventaja de la librería sf, que permite mantener la lógica de data frames al trabajar con objetos espaciales.
Ahora que hemos acotado el área sobre la cuál nos interesa extraer información, utilizaremos la librería osmdata para hacernos de la información necesaria para ejemplificar los diferentes tipos de objetos espaciales que hemos estado estudiando.
El siguiente bloque de código crea el objeto espacial “comercios” mediante la siguiente lógica: en primer lugar se utiliza la función opq() para acceder al servicio de descarga de información de OSM (una API) definiendo el área sobre la cuál se necesitará extraer la información. Esta área la hemos definido con el polígono “area_de_estudio” que representa los límites administrativos de la ciudad de Córdoba. En segundo lugar, utilizamos la función add_osm_feature() para definir qué tipo de datos queremos extraer de OSM. Openstreetmap tiene una manera de catalogar la información que la comunidad carga, definiendo diferentes códigos para diferentes tipos de objetos. Estos códigos pueden tener un primer agrupamiento, dado por el valor “key”, y un segundo agrupamiento dado por el valor “value”. En el ejemplo, nos interesa descargar todos los comercios (shops) en el área de estudio, pero podríamos estar interesados en descarga sólo las panaderías (bakery). En este hipotético caso, deberíamos reemplazar la línea correspondiente por add_osm_feature(key = “shop”, value = “bakery”). Una exploración detallada de los diferentes criterios de catalogación de la información en OSM se puede realizar en el siguiente link.
comercios <- opq(area_de_estudio) %>%
add_osm_feature(key = "shop") %>%
osmdata_sf()Veamos un poco más de cerca el objeto recientemente creado, llamado “comercios”:
comercios## Object of class 'osmdata' with:
## $bbox : -31.5251273,-64.3100961,-31.307545,-64.05725
## $overpass_call : The call submitted to the overpass API
## $meta : metadata including timestamp and version numbers
## $osm_points : 'sf' Simple Features Collection with 3406 points
## $osm_lines : NULL
## $osm_polygons : 'sf' Simple Features Collection with 374 polygons
## $osm_multilines : NULL
## $osm_multipolygons : NULLDe lo anterior se desprende que este objeto contiene 3406 puntos, y 374. Esto se debe a que los usuarios pueden marcar la localización de un comercio con un punto, pero también dibujar un polígono que describa grandes superficies comerciales como los shoppings o ferias. Ya veremos cómo trabajar sobre esta dispersión de criterios, pero de momento, vamos a quedarnos sólo con los 3406 puntos.
comercios = comercios$osm_pointsVeamos un poco en detalle el resultado. La función dim() nos informa sobre la dimensión del data frame resultante.
dim(comercios)## [1] 3406 107Como vemos, el objeto “comercios” tiene 3406 filas y 107 columnas. Es mucha información. Si utilizamos la función names() podemos ver los nombres de cada una de estas columnas:
names(comercios)## [1] "osm_id" "name"
## [3] "addr.city" "addr.country"
## [5] "addr.housenumber" "addr.postcode"
## [7] "addr.street" "air_conditioning"
## [9] "alt_name" "amenity"
## [11] "barrier" "branch"
## [13] "brand" "brand.en"
## [15] "brand.wikidata" "brand.wikipedia"
## [17] "building" "bus"
## [19] "butcher" "check_date"
## [21] "clothes" "contact.email"
## [23] "contact.facebook" "contact.instagram"
## [25] "contact.phone" "contact.website"
## [27] "contact.whatsapp" "currency.ARS"
## [29] "currency.XBT" "delivery"
## [31] "description" "designation"
## [33] "diet.vegetarian" "drink.coffee.togo"
## [35] "drink.wine" "email"
## [37] "entrance" "female"
## [39] "fixme" "healthcare"
## [41] "highway" "instagram"
## [43] "internet_access" "internet_access.fee"
## [45] "is_in" "level"
## [47] "male" "mobile"
## [49] "name.en" "name.es"
## [51] "network" "not.brand.wikidata"
## [53] "note" "official_name"
## [55] "opening_hours" "operator"
## [57] "organic" "outdoor_seating"
## [59] "payment.american_express" "payment.argencard"
## [61] "payment.bank_deposits" "payment.bank_transfer"
## [63] "payment.cabal" "payment.cash"
## [65] "payment.cheque" "payment.contactless"
## [67] "payment.credit_cards" "payment.debit_cards"
## [69] "payment.diners_club" "payment.electronic_purses"
## [71] "payment.gift_card" "payment.italcred"
## [73] "payment.kadicard" "payment.maestro"
## [75] "payment.mastercard" "payment.mastercard_contactless"
## [77] "payment.mastercard_debit" "payment.meal_voucher"
## [79] "payment.naranja" "payment.nativa"
## [81] "payment.provencred" "payment.tarjeta_shopping"
## [83] "payment.visa" "payment.visa_debit"
## [85] "payment.visa_electron" "payment.wire_transfer"
## [87] "phone" "phone.AR"
## [89] "prepayment.cencosud" "public_transport"
## [91] "ref" "ref_code"
## [93] "repair" "second_hand"
## [95] "service" "shop"
## [97] "shop_1" "shop_2"
## [99] "source" "start_date"
## [101] "takeaway" "toilets"
## [103] "trade" "website"
## [105] "wheelchair" "wholesale"
## [107] "geometry"El objeto tiene una columna de geometría (la última) que define las coordenadas geográficas de cada punto, pero también muchas columnas con información adicional, que el usuario que realiza la carga de la información puede o no completar.
Para reducir el ruido de la información resultante, vamos a quedarnos sólo con la columna “shop”, que indica el rubro de negocio al cual se dedica cada uno de los comercios:
comercios = comercios[, c("shop")]
head(comercios)## Simple feature collection with 6 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -64.21934 ymin: -31.43437 xmax: -64.14273 ymax: -31.36709
## Geodetic CRS: WGS 84
## shop geometry
## 643269288 mall POINT (-64.21934 -31.36709)
## 655193217 bakery POINT (-64.18184 -31.42471)
## 663243931 supermarket POINT (-64.1477 -31.42402)
## 664087358 supermarket POINT (-64.14273 -31.43437)
## 705466795 chemist POINT (-64.18341 -31.42552)
## 716005842 alcohol POINT (-64.18664 -31.42056)Y a continuación podemos crear un mapa interativo con los datos resultantes:
mapview(comercios, col.regions = "blue", legend = F)Siguiendo el procedimiento anterior, podemos también extraer otros tipos de objetos espaciales. Por ejemplo, las líneas que representan calles o avenidas catalogadas como “vias secundarias”. Para ello, en la función add_osm_feature() utilizamos la key “highway” y el value “secondary”.
vias_secundarias <- opq(area_de_estudio) %>%
add_osm_feature(key = "highway", value = "secondary") %>%
osmdata_sf()
vias_secundarias = vias_secundarias$osm_lines[, c("name")]
mapview(vias_secundarias, color = "blue", legend = F)## Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, :
## GDAL Error 1: PROJ: proj_as_wkt: DatumEnsemble can only be exported to
## WKT2:2019Por supuesto, también podemos extraer polígonos. En este caso, por ejemplo, extraeremos los polígonos de las escuelas localizadas dentro del área de estudio haciendo uso de la key “amenity” y el value “school”. Entonces:
escuelas <- opq(area_de_estudio) %>%
add_osm_feature(key = "amenity", value = "school") %>%
osmdata_sf()
escuelas = escuelas$osm_polygons[, c("name")]
mapview(escuelas, color = "black", col.regions = "red", legend = F)A continuación, exploraremos algunas estrategias para utilizar datos en formato ráster. Para ello, necesitaremos instalar previamente las librerías terra y geodata.
install.packages("terra")
install.packages("geodata")A continuación, utilizaremos la librería geodata para descargar un modelo digital de elevación en formato ráster. Esta librería permite acceder a informacion libre generada en el marco de la mision Shuttle Radar Topography Mission (SRTM), que determina la elevacion del terreno (en metros sobre el nivel de mar) a traves de un radar incorporado a un satelite. Dado que la informacion se encuentra “precargada” en mozaicos que cubren diferentes extensiones geográficas, se le debe indicar a la función elevation_3s() las coordenadas geográficas de un punto, para que a partir de esta información se puede determinar qué mosaico se va a descargar.
Ahora bien, hasta aquí hemos definido el área de estudio como un polígono definido por el objeto area_de_estudio. Como ya sabemos, este polígono no tiene sólo una, sino varias coordenadas que, leidas en sentido antihorario, definen los límites del área. Y la función requiere de sólo un par de coordenadas, es decir, de un punto. Para sortear este inconvneiente, utilizaremos la función st_centroid() de la librería sf, que permite calcular el centroide de un polígono, entendido como el punto central de un polígono que se calcula a partir del promedio de las coordinadas geográficas que lo definen. Por lo tanto:
coordenadas = st_centroid(area_de_estudio)A continuación, a este nuevo objeto espacial que representa el centroide del polígono del área de estudio, le extraeremos las coordenadas mediante la función de sf st_coordinates:
coordenadas = st_coordinates(coordenadas)
coordenadas## X Y
## [1,] -64.18378 -31.4166En donde X representa la longitud (este-oeste) e Y representa la latitud (norte-sur). Ahora sí, estamos en condiciones de utilizar estas dos coordenadas geográficas para indicarle a la función elevation_3s() de la librería geodata el lugar para el cual queremos descargar la información sobre la elevación del terreno. La función tiene varios arguementos: el primero (lon) requiere que se informe la longitud calculada recién (el valor de X), y para asgnar ese valor tomamos el primer elemento del objeto coordenadas; para informar la latitud (lat) hacemos lo mismo pero con el segundo elemento del objeto coordenadas; finalmente, el argumento “path” requiere informar la ruta en la cual se guardará el archivo a descargar, en donde indicamos que se almacene en un directorio temporal (esto lo pueden modificar ustedes, detallando un directorio específico en sus computadoras):
library(geodata)## Loading required package: terra## terra 1.7.78library(terra)
elevacion <- elevation_3s(lon = coordenadas[[1]], lat = coordenadas[[2]], path=tempdir() )Podemos graficar el ráster descargado:
plot(elevacion)
También podemos hacer un mapa interactivo utilizando la librería mapview, pero aparece un pequeño inconveniente. El objeto “elevación” es de tipo SpatRaster, que es el formato con el que la librería terra trabaja los rásters. Dado que mapview es una librería diseñada con anterioridad a la aparición de la librería terra, la interacción de ambas no es óptima. Para sortear este problema, convertimos el objeto “elevación” de SpatRaster a un ráster clásico de la librería raster.
mapview(raster(elevacion))Sin embargo, el ráster descargado es mucho más amplio que el área de estudio definida. Para adecuar esta información hacemos uso de la función crop() de la librería terra, de la siguiente manera:
cordoba = crop(elevacion, area_de_estudio)
plot(cordoba)
Y al igual que recién, podemos hacer un mapa interactivo con la salvedad de convertir el objeto “cordoba” al formato correpondiente:
mapview(raster(cordoba))Cargar datos de manera local (desde la PC)
Si bien hasta ahora hemos descargado una amplia gama de archivos con datos espaciales directamente desde internet, hay muchas fuentes de información que no tienen sevicios para conectarse directamente y requieren de la descarga previa de algún archivo. Los archivos en formato geográfico más conmunes son aquellos con las extensiones gkpg, kmz, geojson o shp. Sin embargo, cualquier archivo que tenga una columna con la lista de coordenadas geográficas de cada fila puede convertirse en un objeto espacial, incluso una simple hoja de cálculo.
Vamos a utilizar la función st_read() de la librería sf para cargar al entorno de trabajo el archivo barrios_cordoba.gpkg, creando un nuevo objeto llamado “barrios”, y generamos el mapa interactivo correspondiente:
barrios <- st_read("barrios_cordoba.gpkg")## Reading layer `barrios_cordoba' from data source
## `/home/juan/Documentos/Cursos/Ciencia de datos en sociedad/barrios_cordoba.gpkg'
## using driver `GPKG'
## Simple feature collection with 485 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XYZ
## Bounding box: xmin: -64.30993 ymin: -31.52498 xmax: -64.05738 ymax: -31.30852
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84mapview(barrios)Como se mencionó, existen muchos portales que ofrecen descargas de información geográfica que puede ser luego cargada de esta manera al entorno de trabajo en R. Recomendamos revisar los portales poblaciones.org, datos abiertos de la Municipalidad de Córdoba, el portal de la Infraestructura de Datos Espaciales de la Provincia de Córdoba, el portal del INDEC. Y hay muchos otros.
Operaciones con datos espaciales
A continuación, analizaremos algunas funciones de la librería sf que permiten realizar diferentes operaciones con datos espaciales. Comenzaremos por la realización de uniones entre diferentes objetos espaciales. La lógica de esta operación implica, por ejemplo, identificar los atributos de objetos que comparten una localización o localizaciones próximas.Supongamos que nos interesa imputar a cada escuela el barrio al que pertenece. Para ello, utilizaremos la función st_join(). Esta función tiene 3 argumentos. En primer lugar se indica el objeto espacial sobre el que se desea agregar un atributo correspondiente a otro objeto espacial, en este caso, las escuelas. En segundo lugar se indica el objeto espacial desde dónde se desea transferir el atributo correspondiente, en este caso, el nombre del barrio. En tercer lugar, se debe indicar el criterio de union. En este caso, se indica la opción “st_intersects” que establece que ambos objetos deben unirse siempre que se intersecten en algún punto. Hay un cuarto argumento, que no es obligatorio utilizar, pero suele ser útil: se trata de a opción “largest”. En el caso de asignarle el valor “TRUE”, y en el frecuente caso de que una feature de un objeto se intersecte con dos o más features del otro objeto (por ejemplo, que una escuela esté ubicada en el límite de dos barrios, con una parte en uno y otra parte en otro) esta opción vinculará los datos de aquellas que compartan la mayor proporción del espacio.
escuelas = st_join(escuelas, barrios[,c("description")], join = st_intersects, largest = TRUE)## Warning: attribute variables are assumed to be spatially constant throughout
## all geometrieshead(escuelas)## Simple feature collection with 6 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -64.27902 ymin: -31.45692 xmax: -64.1732 ymax: -31.36073
## Geodetic CRS: WGS 84
## name
## 67868340 Instituto Nuestra Señora
## 83552775 Instituto Técnico Salesiano Villada
## 143976290 Escuela Superior de Comercio Manuel Belgrano
## 154249743 Instituto Técnico Salesiano Villada
## 154644719 Colegio San Buenaventura
## 184007376 IPEM 8 Manuel Reyes Reina
## description geometry
## 67868340 AMPLIACION JARDIN ESPINOSA POLYGON ((-64.17479 -31.454...
## 83552775 VALLE ESCONDIDO POLYGON ((-64.27749 -31.363...
## 143976290 ALBERDI POLYGON ((-64.20047 -31.405...
## 154249743 VALLE ESCONDIDO POLYGON ((-64.27666 -31.362...
## 154644719 ALTOS DE SANTA ANA POLYGON ((-64.25112 -31.402...
## 184007376 ALTO VERDE POLYGON ((-64.21076 -31.369...Como vemos, se ha agregado una columna que informa el barrio al que pertenece cada escuela.
Ahora, con esta nueva información, podemos realizar algunas agregaciones no espaciales, pero que tendrán correlato a nivel geográfico. Por ejemplo, supongamos que nos interesa saber la cantidad de escuelas que hay en cada barrio. Para ello, recuperando lo estudiado en módulos anteriores, hacemos uso de las funciones group_by() y summarise() de la librería dplyr. Antes de proceder a calcular ese resumen de cantidad de escuelas por barrio, vamos a remover los atributos espaciales del objeto “escuelas”. Esto lo haremos utilizando la función st_drop_geometry(), que básicamente convierte a data frame todo objeto espacial. Entonces
library(dplyr)
tabla <- st_drop_geometry(escuelas) %>%
group_by(description) %>%
summarise(cantidad = n()) %>%
arrange(desc(cantidad))
head(tabla)## # A tibble: 6 × 2
## description cantidad
## <chr> <int>
## 1 ALBERDI 9
## 2 CERRO DE LAS ROSAS 7
## 3 GENERAL PAZ 7
## 4 ARGUELLO 6
## 5 ARGUELLO NORTE 6
## 6 VILLA AZALAIZ OESTE 6Ahora supongamos que nos interesa confeccionar un mapa que informe sobre la cantidad de escuelas en cada barrio. Vamos a hacer una unión no espacial entre este objeto no espacial llamado “tabla” con el objeto espacial llamado “barrios”. La llave para unir ambas tablas será el campo description que se encuentra presente en ambas. Por lo tanto:
barrios = left_join(barrios, tabla)## Joining with `by = join_by(description)`Ahora sí, confeccionamos un mapa interactivo que informe sobre la cantidad de escuelas en cada barrio de la ciudad. Como es usual, hacemos uso de la función mapview(), agregando ahora el argumento “zcol” que indica sobre qué columna se pretende agregar una paleta de colores a la representación gráfica.
mapview(barrios, zcol = "cantidad")Es importante notar que muchos barrios (331 de un total de 485) tienen NA en el campo “cantidad”. Esto se debe a que, del total de escuelas descargadas desde OSM, ninguna se encuentra localizada dentro del polígono de cada uno de estos barrios. Conviene recordar en este punto que la fuente de información desde la cual se extrajo el dato de escuelas no es oficial, sino que se trata de un mapa colaborativo en el que, seguramente, faltan varias por agregar (aprovachamos para reforzar la necesidad de colaborar con OSM para mejorar la calidad de esta información).
Sobre el mapa anterior podemos agregar las escuelas, para asegurarnos de que, efectivamente, no hay ninguna en los barrios que aparecen con valores NA (en gris en el mapa anterior):
mapview(barrios, zcol = "cantidad") +
mapview(escuelas, col.regions = "red")Otra operación espacial interesante consiste en la unión de diferentes features. Por ejemplo, supongamos que nos interesa unir los barrios Centro, Nueva Córdoba y Alberdi. Lo primero que vamos a hacer es filtrar esos barrios utilizando la función filter().
barrios <- barrios %>% filter(description %in% c("CENTRO","NUEVA CORDOBA","ALBERDI"))
mapview(barrios)Como vemos, el objeto espacial “barrios” tiene ahora sólo 3 filas o features, una para cada barrio. Vamos a unir estas 3 features en un solo polígono, utilizando la función st_union() de la librería sf.
barrios = st_union(barrios)
mapview(barrios)Como vemos, ahora el objeto tiene sólo un polígono, resultante de la union de los 3 barrios anteriores.
Supongamos que nos interesa ver sólo las escuelas que se encuentran dentro de este nuevo polígono. Para ello, utilizaremos la función st_intersection() de la librería sf. Esta función tiene sólo dos argumentos: el primero consiste en la capa sobre la que se quiere seleccionar las features (las escuelas en nuestro ejemplo), y el segundo representa el área sobre la cuál se desea realizar la selección (el polígono de los 3 barrios unidos). Entonces:
escuelas = st_intersection(escuelas, barrios)## Warning: attribute variables are assumed to be spatially constant throughout
## all geometriesmapview(escuelas)El objeto espacial resultante contiene ahora exclusivamente las 15 que se encuentran localizadas dentro de los 3 barrios seleccionados.